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Abstract Escherichia coli is a motile bacterium that moves up a chemoattractant 
gradient by performing a biased random walk composed of alternating runs and 
tumbles. Previous models of run and tumble chemotaxis neglect one or more fea- 
tures of the motion, namely (i) a cell cannot directly detect a chemoattractant gra- 
dient but rather makes temporal comparisons of chemoattractant concentration, 
(ii) rather than being entirely random, tumbles exhibit persistence of direction, 
meaning that the new direction after a tumble is more likely to be in the forward 
hemisphere, and (iii) rotational Brownian motion makes it impossible for an E. 
coli cell to swim in a straight line during a run. This paper presents an analytic 
calculation of the chemotactic drift velocity taking account of (i), (ii) and (iii), 
for weak chemotaxis. The analytic results are verified by Monte Carlo simulation. 
The results reveal a synergy between temporal comparisons and persistence that 
enhances the drift velocity, while rotational Brownian motion reduces the drift 
velocity. 
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1 Introduction 

As described by Berg [Q], Escherichia coli is a common intestinal bacterium with 
a body m 1 /im in diameter and 2 /im long. A typical E. coli cell possesses 
approximately 6 helical fiagella, each w 20 nm thick and « 10 fim long, which 
emerge from random points on the cell membrane. Each fiagellum is powered by 
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a reversible rotary motor. When the flagella spin counter-clockwise, they form a 
synchronous bundle and propel the cell in an approximately straight path called a 
'run', with a swimming velocity of 20 jiims" 1 and a corresponding Reynolds 
number of Re w 10~ 5 . The run is not entirely straight, as the cell is subject to 
rotational and translational Brownian motion due to thermal collisions with the 
molecules in the surrounding fluid. Translational Brownian motion is negligible 
compared to the cell's swimming motion, but rotational Brownian motion is sig- 
nificant and in water at room temperature it causes a root-mean-square angular 
deviation in radians of about 0.5 f 1 / 2 , where t is in seconds. When one or more 
flagella spin clockwise, the bundle comes apart and the cell moves in an erratic 
path with little net displacement called a 'tumble', which reorients the cell to face 
a random new direction. Tumbles exhibit persistence of direction (hereafter sim- 
ply referred to as 'persistence'), meaning that the direction faced by the cell after 
a tumble is not isotropically random. Rather, the mean angle between the new and 
the previous direction is 62° [Oil, so the new direction after a tumble is more likely 
to be in the forward hemisphere. 

An E. coli cell performs chemotaxis by executing a biased random walk com- 
posed of alternating runs and tumbles. The run durations are well described by 
a Poisson interval distribution (i.e. exponential distribution) with a rate constant 
Qumble rate') of A w 1 s -1 , and a corresponding mean run duration of A -1 ~ 1 s 
yfl. Tumble durations also follow a Poisson interval distribution, but are much 
shorter than runs, with a mean duration of « 0.1 s. An E. coli cell is too small to 
detect spatial differences in the concentration of a chemoattractant on the scale of 
a cell length. Instead, the cell performs temporal comparisons. As the cell swims 
about it continually measures the concentration of chemoattractants in its environ- 
ment (e.g. serine, aspartate) and compares the chemoattractant concentration over 
the past second with the concentration over the previous three seconds. The cell 
modulates the tumble rate in response to a concentration difference, so that runs 
up the chemoattractant gradient are extended. The direction-dependent tumble rate 
causes the cell to drift toward regions of high chemoattractant concentration. 

One measure of the effectiveness of chemotaxis is the drift velocity, defined 
as the mean velocity at which a cell moves up a chemoattractant gradient. No pre- 
viously published model of chemotaxis allows the calculation of the chemotactic 
drift velocity while simultaneously accounting for temporal comparisons, persis- 
tence, and rotational Brownian motion. Existing literature on the topic includes 
the following. Patlak [ 10] presented an analysis of random walks with persistence 
and bias, but his framework does not accommodate temporal comparisons. Brown 
and Berg [5] performed Monte Carlo simulations for a model in which a cell pos- 
sesses an exponentially decaying memory; they included persistence but not rota- 
tional Brownian motion. Schnitzer lfl2ll presented a theory of random walks with 
persistence and rotational Brownian motion, but without temporal comparisons, 
de Gennes [7] presented an analytic method for calculating the chemotactic drift 
velocity taking account of temporal comparisons, but neglected Brownian motion 
and persistence. Recently, Erban and Othmer [8|] achieved the noteworthy task of 
deriving a complete advection-diffusion equation for a suspension of chemotactic 
organisms taking account of both temporal comparisons and persistence. 
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In this paper we extend the model of de Gennes [7] to include Brownian mo- 
tion and persistence, and we report a new effect: persistence can markedly enhance 
the drift velocity. 

2 The model 

Consider a cell performing run and tumble chemotaxis, swimming in an unbounded, 
stationary fluid which contains a uniform chemoattractant concentration gradient 
Vc || z, where z denotes the unit vector in the z direction. During a run, the cells 
swims at constant speed v s , and the probability that the cell tumbles in the next 
time interval dt is X(t)dt, where X is the 'tumble rate'. Like de Gennes I2D, we 
assume that for weak concentration gradients the tumble rate depends linearly on 
chemoattractant concentration history, so that the tumble rate is given by: 

X(t)=Xo[l-A(t)}, (1) 

where Xq = Is -1 is the baseline tumble rate, and the fractional change in tumble 
rate is given by 

A{t)= J c{t')R(t-t')dt', (2) 

where c(f') is the chemoattractant concentration experienced by the cell at time 
t', and R is the cell's 'response function'. Note that R may be thought of as the 
impulse response of the tumble rate, since it describes the way that X changes 
when the cell is subject to a Dirac delta impulse of chemoattractant concentration. 
Equation (Q~|) also appears in the papers of Schnitzer lfl2ll and Clark and Grant 
lHO, and it is motivated by the experimental results of Block et al 0] and Segall 
et al 10]. Our analysis will be restricted to 'weak chemotaxis', meaning small 
fractional changes in the tumble rate, i.e. \A (t) | <C 1 . Physically, weak chemotaxis 
corresponds to a shallow chemoattractant gradient. 

Rotational Brownian motion causes the cell to gradually veer off course, so 
that in between tumbles the probability density function / of the swimming direc- 
tion e evolves according to the Fokker-Planck equation 

it=D R V 2 J, (3) 

where D R is the rotational diffusion coefficient and is the Laplacian in direction 
space. Berg [2] estimated D R w 0.062 s^ 1 for an E. coli cell swimming in water 
at room temperature. 

When a cell tumbles, its choice of new direction is governed by a probability 
distribution which is axisymmetric about the initial direction. We allow for the 
tumbles to exhibit directional persistence, so that the expected scalar product of 
the swimming directions e(0~) and e(0 + ) immediately before and after a tumble 
at time f = is given by 

E[e(0-)-e(0 + )] = a p , (4) 

where Op is the 'persistence parameter'. Experimentally, Op « 0.33 J2]. 

The drift velocity, v<i, can be defined in terms of the expected motion in a single 
run. Let Zf be the z location of a cell at the end of a run, relative to its position at 
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the beginning of a run, and let tf be the duration of a run. We treat the duration of 
tumbles as negligible, so 

v d = E[ Zf ]/E[t f ] (5) 

where E denotes an expectation value. If the non-zero duration of tumbles is taken 
into account, then the only effect is to reduce Vd by a small factor. In taking expec- 
tation values, we assume that the cell has already been swimming in the fluid 
for a sufficiently long time that the probability density function of Zf is time- 
independent. While our model assumes an unbounded domain, we note that for 
a cell swimming in a bounded domain of length 3> v s /Ao, provides a measure 
of the transient average velocity up the chemoattractant gradient before the cell 
encounters the boundaries. 



3 Calculation of drift velocity for general response function 

Calculating is non-trivial because of the interdependence of the tumble rate and 
the path taken by the cell. The tumble rate at any time depends in principle on the 
entire path history of the cell through (0, while the path of the cell depends in 
turn on the tumble rate. As noted above, in order to make the analysis tractable, 
our analysis will be restricted to weak chemotaxis, i.e. < \A (t) | -C 1. In this case, 
one has E[z f ] = v s Ao<9(A) and E[t f ] = [1 + 0(A)}/^, so 

v d = Ao£[z f ]+v s O(zi 2 ). (6) 

We shall neglect terms that are 0(A 2 ). 

Consider a run commencing at time f = at location z = 0. During the run, the 
cell swims in a random walk governed by rotational Brownian motion until the 
run terminates with a tumble at time tf. Note that the expected stopping location 
of a terminated random walk with a stopping rate A sto p is the same as the expected 
first event location on an unterminated walk with an event rate /l event = A sto p . Thus, 
in calculating E[zf], it is permissible to treat tumbles for f > as events that have 
no effect on the cell's motion, and treat the tumble at t = tf, z = Zf as a first event 
(the first tumble in t > 0). The utility of this treatment is that we may conceptually 
break the expectation E in E[zf] into two consecutive operations. First, assuming 
a given path z(t) : — °° <t<°° taken by the cell, one calculates the conditional 
expectation of Zf for that path. Second, one takes the expectation over all such 
paths to obtain E [zf] , with the understanding that in the t < section of a path the 
cell is subject to reorientations due to both Brownian motion and tumbles, whereas 
in the t > section of a path the cell is subject to reorientation due to Brownian 
motion alone. Writing out the two expectations in symbolic notation, 



E [Zf] — -Epaths 



dtZ p - dt h(t)Ppath(t) 



(7) 



where the £ pa ths denotes an expectation over paths, z pa th(0 denotes the position 
of the cell at time I on a particular path, and /? pa th(0 1S the probability density 
function for the tumble time tf on a particular path. Since tumbles for / > are 
treated as having no effect on cell motion, paths are independent of tf and one is 
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free to take the path expectation inside the integral over tumble times. Dropping 
the 'path' subscript for brevity, one then has 



E[z f ]= / dtE[z(t)p(t)], (8) 
Jo 

The probability density function p(t) for the tumble time is given by 



p(0=A(f)exp 



- / X(t')dt' 



JO 



(9) 



where /1(f) is the path-dependent tumble rate at time f. Substituting © into (|8]l 
and © and integrating by parts yields 



Vd 



Jo 



dtE 



->(f)exp 



Mt')dt' 



(10) 



where w(t) =dz(t) /dt. Writing X{t) = Ao[l — 4(f)] in ( flOl) . expanding the expo- 
nential in powers of A and keeping only the linear term, one obtains 



Vd : 



Jo 



dtE 



w(t)e 



1+Ao [ A(t')dt' 
Jo 



(11) 



We shall calculate the drift velocity for the case where R is given by a Dirac 
delta function and later generalise to an arbitrary response function. Consider the 
response function 

R(t)=AS(t-T). (12) 
With this R, the fractional change in tumble rate is simply 



A{t) =Ac(t-T)} 

= Ac Q +A\Vc\z(t-T)}, 



(13) 



where c(t) =cq + \ Vc|z(f ) is the chemoattractant concentration seen by the cell at 
time t. We shall discard the Co term, since eventually we shall require the response 
function R to have zero mean, so that an additive constant to c in (0 has no effect 
on A. Substituting ([T3T l into dTTb and discarding terms containing cq, one finds 



v a = Xq( dte-*E\w{t)]+l%A\Vc\( dteT^ [ dt' E[w(t)z(t' -T)]. (14) 
Jo Jo Jo 

We may write 

E[w(t)]=v s z-E[e(t)]. (15) 

The probability density function for the angular deviation caused by rotational 
Brownian motion is symmetric about initial swimming direction, so we may write 

E[e(t)]=E[e(t)-e(0 + )]E[e(0 + )} (16) 

where e(0 + ) is the cell's swimming direction at the beginning of the run. One can 
straightforwardly solve the Fokker-Planck equation |3]) to find that the direction 
correlation function in the absence of tumbles is 

E[e( tl )-e(t2)}=c- 2DR \'^ (17) 
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(t\,t2 > 0), and substitution of dT71) into dT6b and (f]~5T > gives 

E[w{t)]=c' 2DR, E[w(0+)] 1 (18) 

where w(0 + ) is the z component of the cell's velocity at the beginning of the run. 
Substituting ( 118b into (114b and swapping the order of integration over f and t' 
yields 

v d = £[w(0 + )]Ao/(Ao + 2Z) R ) 

+XoA\Vc\ f dt' f dte~^E[w(t)z(t'-T)] (19) 
Jo Jt' 

To proceed further, we must find expressions for £[w(0 + )] and £[w(*b)z(*a)]> 
for — oo < ? a < ^ and < t b < °°. Since we are neglecting terms 0(A 2 ) in vj, 
and the A\ Vc| factor in ( 119b comes from one power of zi, we can neglect chemo- 
taxis altogether in calculating E[w(t\,)z(t a )]. Neglecting chemotaxis, position and 
velocity are governed by isotropic distributions, so 



E[w(t b )z(t a )] = / dtE[w(t b )w(t)} 
Jo 

= \vl [* dtE[e(t h )-e(t)}. (20) 
i Jo 

For < ? a < t\,, we can use equation dl7l ) for the direction correlation function to 
find that 

E[ W {t h )z{h)] = \vl tdte 2DR ^ 
j Jo 

= ^e- 2DRtb (e 2DR '°-l),0<t a <t b . (21) 

For f a < [and hence t < in d20ll, we expand the direction correlation function 
into a product of direction correlations between different times, and then simplify 
using © and ( fTTt : 

£[«(fb)-«(f)] = «[«(«b) ^(0 + )]£[e(0 + ) .«(0-)]£[«(0") •«(*)] 

= cepe- 2DR ^£[<?(0-)-<?(0]. (22) 

To obtain an expression for Zs[e(0~) we note that for? < the cell undergoes 
reorientations due to both rotational Brownian motion and tumbles. Lovely and 
Dahlquist [9] showed that the direction correlation function for isotropic run and 
tumble random walks without rotational Brownian motion is 

E[e{t x )-e{t 2 )}=^-^\ (23) 

where 

1 (24) 



Ao(l-Op) 

is the direction correlation time. Lovely and Dahlquist |@] also showed that rota- 
tional Brownian motion can be regarded as a Poisson tumbling process, and that 
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if two Poisson tumbling processes have individual direction correlation times X\ 
and Vi, then the direction correlation for the combined process is still of the form 
(123 b . but with correlation time T c given by 

I-I + I. (25) 

T c Ti T 2 

One may think of T c as the characteristic time it takes the cell to perform a 'U- 
turn' . We apply d24l l and d25T l to our problem to find that 

T c = l/[Ao(l-Op)+2Z) R ] (26) 

and 

E[e(0-)-e{t)}=c^- a ^ +2D ^. (27) 
Substitution of d27j and $22^ into (|20| yields 

EWMti] = ^ S 3^(1-^+2^] 1 > ^ * (28) 

Having found E[w(t b )z(t a )], we turn our attention to £[w(0 + )]. In the absence 
of rotational Brownian motion, Z?[w(0 + )] = 0, since although the cell alters run 
durations in response to the chemoattractant, run directions remain random and 
isotropic. In the presence of Brownian motion, however, this is no longer true. 
During a run, rotational Brownian motion causes the cell's swimming direction 
to drift. For instance, modulation of the tumble rate might make it more likely 
that the cell tumbles when heading down the chemoattractant gradient, so there 
is down-gradient bias in the expected swimming direction at the end of a run. If 
0£p > the cell is consequently more likely to commence each run heading down- 
gradient. 

From the definition (|4| of the persistence parameter, we may write 

E[w(0 + )] = a p E[w(0-)}. (29) 

Furthermore, if we assume that the cell has been swimming in the chemoattrac- 
tant gradient for a time f sw j m 3> T c [where x c is given by J26H . then the expected 
swimming velocity just before a tumble does not change from one tumble to the 
next, so £[w(ff - )] = i?[w(0~)] and hence 

E[ w (0 + )] = a p E[w(t f -)}. (30) 

As in the derivation of ©, we may break the expectation on the right hand side of 
(130b into an integral over tumble times and an expectation over paths: 

E[w(tf)}= f dtE[w{t)p{t)} (31) 
Jo 

Using (O and (fl]) to expand p(t) in powers of A and keeping only the linear term, 
one obtains 



p(t)=Xoe 



1 



-A|Vc|z(f-r) + AoA|Vc| f z(t'-T)dt' 
Jo 



(32) 
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Substituting d32l into d3ll and re-arranging yields 

E[w(tf)]=X l~ dt€-*E[w{t)]-koA\Vc\{h-h), (33) 
Jo 



where 



and 



h= [ dte-**'E[w(t),z(t-T)], (34) 
Jo 



/ 2 = Ao (~ dt' rdte-^E[w(t),z(t'-T)]. (35) 
Jo Jt' 

Using ( fT8l ). 031 ) simplifies to 

E[w{tf)\=E[w{0 + )]X /{^ + 2D R )-X A\yc\{h-h). (36) 
Eliminating £[w(ff~)] from d3Qb and d36b then yields 

£W0+)] = w»^). 

Using the definitions of I\ and 7 2 and equations ( 1281) and ( 1211) for £[w(fb),z(f a )]> 
one finally obtains the expression for Z?[w(0 + )]: 

2A|Vc|a p D R v^[(A <) + 2D R )e-^( 1 -»p)+ 2oR ] r -Aoa p ] 
£[W(0 )] = 3(A + 2D R )[Ao(l-a p )+2Z) R ]2 ' (38) 

We now have all the necessary pieces to write down the expression for Vd- 
Substituting < l38T l, d28l l andtETTl into ( TT9l l, we obtain: 

v d =^(r), (39) 

where 

v 2 |Vc|A e-( A o +2Z)R ) r (l-ap)[(Ao + 2D R )e^ 0i P r -A{ ) ape^ +2oR ) r ] 
( } ~ 3(Ao + 2Z)R)[A<,(l-ap)+2DR]2 

(40) 

The function is effectively the Greens function for an impulse response function. 
For a general response function 7?, the drift velocity is 

v d = r°dTR(T)k(T). (41) 
Jo 

Our calculations have been based on the assumption of a uniform concen- 
tration gradient Vc. In fact, our results still apply to the case of a non-uniform 
gradient, provided that the distance over which Vc varies is much larger than the 
distance covered by the cell in a time T c , i.e provided that 

^ »v s t c . (42) 



Similarly, our calculations are based on the assumption that c has no time de- 
pendence in the Eulerian reference frame. However, our results are still valid for 
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time-dependent c, provided that the rate of change of c seen by a stationary cell is 
much slower than the rate of change of c seen by a moving cell i.e. provided that 



dc 
dt 



< v s |Vc| 



(43) 



Calculation of Vd for time-dependent c might be relevant, for instance, to problems 
where the chemoattractant field is altered by cells themselves through secretion or 
degradation. 



4 Consistency with known results 

de Gennes [7] derived an expression for the drift velocity in the absence of rota- 
tional Brownian motion and persistence. Setting D R = and Op = in ( 140b and 
(14 X b one finds 



v 

Vd = 



2 'Vc| 



[ dTR(T)e~^ J ', (44) 
Jo 



3Ao 

consistent with equation 16 of [Q]. 

Schnitzer derived expressions for the flux of cells in a suspension of non- 
interacting chemotactic cells under the assumption that each cell modifies its tum- 
ble rate according to its instantaneous swimming direction relative to the chemoat- 
tractant gradient, so that 

X = l-eee z . (45) 
We note that within our framework ( |43T > is equivalent to setting 

|Vc|v s AT 

and taking the simultaneous limits A T — > and T — > 0. Substituting (l46l l into d4Tb 
and taking A T — > and T — > yields 

_ e^o(l-qp) 
Vd 3[2 J DR + Ao(l-a p )]' l4/j 

The flux J of cells in a uniform suspension is simply proportional to the drift 
velocity. For cells swimming with unit velocity, 

eAo(l-Op) 

/= 3[2^ + Ml-ap)] Pg " (48) 

where p is the number density of cells. Equation (148b is consistent with equations 
6.8 (dp >0,D R = 0) and 7.6 {D R > 0, Op = 0)of I12ll in the case of a uniform 
suspension. In passing, we note that Schnitzer 11211 did not consider the case where 
both D R and Op are non-zero; in this case d48l) shows that if D R > then |/| is a 
decreasing function of Op and |/| — > as Op — > 1. 

The work of Erban and Othmer |§] contains a derivation of an advection dif- 
fusion equation for a suspension of chemotactic organisms, allowing for temporal 
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comparisons and persistence. Temporal comparisons are modelled using internal 
states of the cell y\ and yi that evolve according to 

. f(c(t))-yi-y 2 

• Vl = k 

n = (49) 

'a 

where t e is an 'excitation time' and f a > t e is an 'adaptation time', and f(c) is the 
fraction of the cell's receptors bound to chemoattractant molecules. The tumble 
rate is set as A = Xq — by\ . 

The Erban and Othmer [8] model for A (f) can be cast in the same form as ([T|) 
and (0 by finding the general solution to ( 149b : 

yi(0 = f_Jt'f(e{t'))Y{t-t'), (50) 

where 

y( " = ^(— J (51) 

is the Green's function for yi. For small changes in concentration around Co, one 
then finds that the Erban and Othmer JH] model for X (t) is equivalent to our model 
with a response function 

R(t) = ^f'(c )Y(t). (52) 
Substitution of (l52l into (HTT l then yields a drift velocity of 

= bt.y s f(c Q )\vc\ 

d 31o[l + (l-a p )Va][l + (l-Op)Ve]' ^ ' 

which is identical to the drift velocity that one can calculate from equations 4.37 
and 4.40 of H. 



5 Drift velocity for a 'realistic' response function 

5.1 Choice of response function 

The validity of © has never been directly assessed, and the chemotactic response 
function R has never been directly measured. However, Block et al [4] and Segall 
et al [ 14] monitored the response of a single rotary motor on an E. coli cell by teth- 
ering the flagellum to a glass surface and delivering small impulses of chemoat- 
tractants to the cell's immediate environment. By repeating the experiment mul- 
tiple times and with different cells, the experimenters made measurements of the 
motor bias (the probability of counter-clockwise rotation) as a function of time, 
(see for example figure 1 of Segall et al ifbUl ). The impulse response of the mo- 
tor bias is double-lobed; the bias is raised above the baseline for the first wis 
after the delivery of the impulse, reduced below the baseline for the following 
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Fig. 1 The response function used in the results section of this paper, originally derived by Clark 
and Grant QJ. 



« 3 s, and then returns to baseline. The areas of the two lobes of the response 
are equal. Furthermore, the experimenters found that the responses to other time- 
series of stimuli (e.g. ramp or sinusoid changes in chemoattractant concentration) 
are consistent with the cell behaving as a linear system, so that the response to 
an arbitrary stimulus is well described by the convolution integral of the stimulus 
with the impulse response (c.f. equation[T|. The primary exception to the linear be- 
haviour is that for small changes in chemoattractant concentrations, cells respond 
to increases in concentration but not decreases; we neglect this nonlinearity in our 
analysis. 

In the absence of experimental measurements of R, a convenient assumption 
is that R has a similar shape to the impulse response of the individual motor bias 
reported in Block et al Q40 and Segall et al 11411 . scaled and shifted appropriately 
[121. This means that (i) R should be composed of a positive lobe followed by a 
negative lobe of equal area so that 

[ dtR(t)=0, (54) 
Jo 

and (ii) R(t ) should decay to zero for t greater than about 4 s. (Note that the equiv- 
alent response function ( 1521) of the Erban and Othmer [8] model always satisfies 
(i) and also satisfies (ii) provided that f a < 4s.) Additionally, we require \A (t) | <C 1 
in order for our linear analysis to be valid. For R satisfying ( 154b . \A(t) \ is largest 
when the cell swims straight up or down the chemoattractant gradient, with the 
maximum possible value of \A (t ) | given by 

|4(0|max = V s |Vc| / dttR{t). (55) 

Jo 

We choose R to be 

Of! 2 

R(t) = e -«ot[l _ Aot/2- (Aot/2) 2 ], (56) 

for which |4(f)| max = £, with £ -C 1 in order for our linear analysis to be valid. 
This R, plotted in figureQ] is a theoretically motivated response function derived in 
|0] that has the above-mentioned properties (i) and (ii), and it matches the experi- 
mentally measured motor bias reasonably well. It may seem odd that |Vc| appears 
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Fig. 2 Drift velocity vj as a function of persistence Dtp, for rotational diffusivities D R = 
0,0.01,0.062,0.2 radians 2 s~' (top to bottom), with Ao = Is -1 , if- = 0.062 radians 2 s~' is an 
estimate for an E. coli cell swimming in water at room temperature |2|]. 
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Fig. 3 Drift velocity vj as a function of persistence Op, for rotational diffusivity D = 
0.062 radians 2 s~' with Aq = ls~'. The dashed curve is for a cell that responds to the instan- 
taneous chemoattractant gradient [equation d47H and the solid curve is for a cell that performs 
temporal comparisons of chemoattractant concentration requationj57H. 



on the right hand side of (|56l ). since R is a property of the cell rather than its envi- 
ronment. However the factor of 1/| Vc| in R simply reflects the fact that we have 
used £ to parameterise the strength of the combined effect of the cell's response R 
and the chemoattractant gradient Vc on A (t ) . 



5.2 Dependence of drift velocity on parameters 

Substitution of (|56l l into (HTJ yields 

£v s A 3 [Ao(5 -2ap) +4£> R ](1 - Op) 
Vd 9[2 J DR + A (l-a p )][2 J DR + A (2-a p )] 3 ' P) 

Figure [2] shows Vd as a function of Op for a variety of values of D R . 

A striking feature of figure [2] is that the maximum of Vd occurs at a positive 
value of Op, so positive persistence increases Vd (up to a point). This is markedly 
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different to the case of a cell responding to the instantaneous chemoattractant con- 
centration gradient (the case analysed by Schnitzer LI 2]), for which Vd is highest 
for a p = —1. Figure [3] compares the two cases. A qualitative explanation of the 
enhancement of Vd by persistence for a cell performing temporal comparisons is 
as follows. Inspection of (0 reveals that the largest modulations in A occur when 
the cell swims straight up or down the concentration gradient for 4 s or more. For 
larger Op, the cell is more likely to swim approximately straight up or down the 
gradient for 4 s or more (since it undergoes a smaller change of direction when 
it tumbles), hence A is on average modulated by a larger amount, and on average 
the chemotactic response is stronger. Alternatively, one can explain the effect in 
terms of the relevance of past concentration information. If Op is close to 1, then 
information about the concentration over the last 4 s, on the basis of which the cell 
biases its tumble rate, is more relevant than if (Xp = 0. 

Figure also shows that for D R > 0, v d — > as Op — > 1 . This makes sense 
intuitively since for 0Cp — 1 the cell effectively does not tumble at all and hence is 
unable to bias its random walk. However, for D R = it appears that Vd is largest 
when a p = 1, whereas sensibly one expects that Vd = when (Xp = 1. This ap- 
parent paradox can be resolved if one considers the time-scales over which our 
results are valid. In deriving (137b , we assumed that the probability distribution 
governing £[wj~] is steady state. This is true only if the cell has been swimming in 
the chemoattractant gradient for a time f sw ; m ;§> T c , where the direction correlation 
time t c is given by d26l . For f sw j m 3> T c , tumbles and rotational Brownian mo- 
tion have randomised the cell's swimming direction so that the current swimming 
direction is independent of any initial conditions. Clearly T c — > °° as (Xp — > 1 for 
D R = 0, so (Xp — 1, D R = is singularity at which our results are not valid. 

The observed persistence in wild-type E. coli is Op ~ 0.33 and d57l l pre- 
dicts that for a rotational diffusivity of D R = 0.062 radians 2 s _1 , Vd is 38% larger 
for (Xp = 0.33 than for (Xp = 0. Thus, wild-type E. coli cells enhance their drift 
velocity considerably by having non-zero persistence. However, the drift velocity 
is maximised by a persistence of (Xp « 0.78, which is considerably larger than the 
observed persistence of (Xp ~ 0.33. This discrepancy seems surprising, since one 
might expect evolution to have optimised a p to maximise Vd- One possible ex- 
planation for the discrepancy is a compromise between transient and steady-state 
behaviours of the cell. Clark and Grant [6] suggested that in order for a bacterium 
to be successful it must both both (a) move toward nutrients in the short run (max- 
imise Vd) and (b) in the long run spend most of its time near the regions of highest 
nutrient density if confined to a bounded domain. Clark and Grant [6] argued that 
the double lobed form of R can be understood in terms of a compromise between 
(a) and (b), with (a) favouring an R with only a positive lobe and (b) favouring an 
R with only a negative lobe. A similar compromise may govern (Xp. 

Steady state performance is optimised by an (Xp close to zero, for the following 
reason. The translational self-diffusion coefficient of a cell performing runs and 
tumbles is D = t c v^ /3 19] and, from ( 1261 ), D is a monotonically increasing function 
of (Xp. Now, suppose a cell is confined to a bounded domain with uniform Vc. 
Ideally, in the long term, the cell would spend all its time right at the boundary 
with the highest c. In fact, the cell wanders within a characteristic length L of that 
boundary, where L w D/vd is governed by a balance between diffusion away from 
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the boundary and drift velocity toward the boundary. For D R = 0.062 radians 2 s _1 , 
one can show that L is minimised by CCp ~ 0.06. A similar argument holds if one 
considers a distribution of chemoattractant with a local maximum, with the same 
conclusion that long term performance is optimised by Op close to zero. Thus, the 
observed persistence of Op « 0.33 might reflect a compromise between transient 
and steady-state performance. 

A second factor that disfavours large Op is that E. coli are denser than water, 
which may cause them to eventually swim downward in the absence of strong 
tumbles. For instance, Salmonella typhimurium is a denser-than-water bacterium 
that performs run and tumble chemotaxis in a similar manner to E. coli, and mutant 
non-tumbling S. typhimurium cells swim downward on average and collect at the 
bottom of the vessel in which they are contained [1]. This is presumably because 
the hydrodynamic drag on the flagella is larger than that on the cell body, so a 
torque is exerted on a non-vertically-oriented cell Hill . The implication is that it 
would be disadvantageous for a cell to have Op very close to 1 since it would then 
swim downward on average, regardless of the direction of the chemoattractant 
gradient. 

Other factors potentially influencing the optimal value of Op for a real E. coli 
cell are that (i) an environment with background fluid flows or with rapidly chang- 
ing distributions of nutrients may favour a smaller <Xp, so that the cell may perform 
a 'U-turn' more quickly, and that (ii) non-linear effects have been neglected in our 
analysis and may play an important role for E. coli in its native environment. 



6 Confirmation of results by Monte Carlo simulation 

Monte Carlo simulation was used to verify the correctness of the analytic calcula- 
tions for £ 1 and to assess the accuracy of the analytic calculations for larger £. 
For the purpose of simulations, we used a different response function: 



This form of R was chosen for computational convenience; only the most recent 
four seconds of a cell's trajectory need to be stored in memory. Like the R defined 
by ( |561 l, the R defined by d58l l possesses a positive lobe followed by a negative 
lobe of equal area, and \A (f ) | max = £. Note that we do not have to specify a value 
for | Vc|, since R and | Vc| appear together as a product in the equation for A (t). 

The simulation method was as follows. For each set of parameter values, 100 
simulations were run, each simulating 10 5 s of swimming. Time was discretised 
into steps of size At = 0.01 s. At each time-step, the tumble probability was cal- 
culated using a discrete approximation of the integral (O, and a random number 
generator was used to decide whether the cell tumbled. If the cell tumbled, then the 
new direction faced after the tumble was chosen from an axisymmetric distribution 
about the old direction, such that the new direction made an angle of arccos(Op) 
with the old direction. If the cell didn't tumble, then rotational Brownian motion 
was simulated by giving the cell a new direction chosen from an axisymmetric 
distribution about the old direction, such that the new direction made an angle of 




(58) 
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Fig. 4 Drift velocity as a function of persistence a p for Z) R = (top) and Z) R = 
0.062radians 2 s~ 1 (bottom), with e = 0. 1 in both cases. Solid lines show analytic results, squares 
show simulation data, and error bars show 95% confidence intervals. 
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Fig. 5 Drift velocity vj as a function of e with D R = 0.062 radians 2 s~' and a p = 0.33. The 
solid line shows the analytic prediction and squares show simulation data. The 95% confidence 
intervals are smaller than the squares. Note that the analytic predictions are within 10% of the 
simulation data for e < 0.6. 



arccos(l — 2D R At) &s 2\/D R At with the old direction (note that this is consistent 
with equation W7\. For each simulation, the drift velocity was estimated as the net 
Z displacement divided by the simulation duration. The baseline tumble rate was 
set to Ao = Is -1 for all simulations. 

For e = 0. 1, there is excellent agreement between the simulation data and an- 
alytic calculations, as shown in figure [4] For larger £, the accuracy of the analytic 
calculations is diminished. Figure [5] shows that the results of the analytic calcula- 
tions are within 10% of the simulation results for £ up to approximately 0.6. 

7 Discussion and Conclusion 

In this paper we presented an analytic method for calculating the drift velocity of 
a cell performing run and tumble chemotaxis, taking account of (i) temporal com- 
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parisons, (ii) persistence, and (iii) rotational Brownian motion, and we verified 
the results with Monte Carlo simulations. The calculations are novel in that they 
are the first calculation of chemotactic drift velocity simultaneously to include all 
features (i), (ii) and (iii) of the motion. Using estimates for the chemotactic re- 
sponse function and the rotational diffusivity, our key finding is that persistence 
can markedly increase the drift velocity. For instance, a persistence of CCp = 0.33 
(characteristic of wild-type E. coli |2j]) can increase the drift velocity by « 40% 
relative to Op = 0, while a persistence of Op = 0.78 can double the drift velocity 
relative to Op = 0. 

There are a number of possible reasons for discrepancy between the value of 
Op that maximises and the value of Op observed in wild-type E. coli. First, 
our model looks only at behaviour in an unbounded domain, or equivalently at 
transient behaviour in a bounded domain, whereas one expects that in order to 
be successful a bacterium must optimise both its transient and its steady state be- 
haviour. Indeed, it has been proposed that the double-lobed form of R reflects a 
compromise between transient and steady state behaviour [6]. The same might be 
true of Op, with a smaller value of Op resulting in a more favourable steady-state 
behaviour. Second, our analysis neglects the fact that E. coli cells are denser than 
the fluid in which they swim and geometrically asymmetric (front-to-back), so that 
a cell with Op ~ 1 would swim downward on average. Third, our analysis does not 
include rapidly varying chemoattractant distributions or background fluid flows. 
Finally, our analysis is linear and therefore limited to weak chemoattractant gra- 
dients; it is not inconceivable that the dependence of drift velocity on persistence 
might be different in stronger chemoattractant gradients. 
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A Solving the Fokker-Planck equation 

The Fokker-Planck equation is 

J-=D R V 2 /, (59) 

where is the Laplacian in direction space and e is the swimming direction. We can describe 
e by spherical polar coordinates, so that e ■ e z = cos 6 and f(6,tj>) sm6d6d(j> is the probability 
that e lies in the solid angle sin Q d& dtj> . Assume axisymmetry, so that / is a function of 8 only. 
Writing out the Laplacian in full, we then have 

^WS + cote^Y (60) 



dt \de 2 de 

Without loss of generality, we choose our initial condition to be e(0) = e-. Since we are only 
interested in the direction correlation function, we need only calculate £[cos0(?)] rather than 
solve for /. Consider 

f df 

dE[cos6]/dt = 2n sing cos0 -±d0. (61) 
Jo at 

Substituting 460b into j6Tb and integrating by parts, one finds 

dE[cos6]/dt = -2D R {cos6), (62) 

so 

£[cose(0]=e- 2flR '. (63) 
Generalising from this result it is straightforward to show that 

E[e(t l )-e(t 2 )]=e- 2DR ^. (64) 



